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1 Introduction 

Coupled second order ordinary linear differential equations, coupled oscillators for short, play 
an important role in almost all areas of science and technology (see the introduction of [1] for a 
recent review). The phenomena of coupled systems appear on all length- and time-scales: from 
synchronization of power generators in power-grid networks [2,3], through the traffic control 
of vehicular platoons [4,5,6,7,8,9], collective decision-making in biological systems [10,11,12, 
13,14] (e.g., transfer of long-range information in flocks of birds), to the atomic scale lattice 
vibrations (so-called phonons), just to name few of them. The nature of communication within 
such a systems crucially influences the behavior of it. In the presence of centralized information, 
e.g., the knowledge of the desired velocity by members of a flock, the performance of many of 
these systems is good [4,7,8] in the sense that the trajectories of the agents quickly converge to 
coherent (or synchronized) motion. On the other hand, in decentralized systems convergence to 
coherent motion is much less obvious, since no overall goal is observed by all agents. In this case, 
the only available observations (i.e., of position and/or velocity) are relative to the agent. The 
complication of the problem is even greater if information is exchanged only locally - by agents 
in a neighborhood that is small in comparison to the system size. 

As an aside, we point out that there exist another class of somewhat similar problems, also 
with a wide range of applications, namely the dynamics of consensus^ see [15,16]. The difference 
is that in our case the agents are Newtonian (are subject to force mx, i.e., mass x accelerations), 
whereas in consensus type equation, that is not the case. Therefore consensus equations tend to 
be coupled first order ordinary differential equations with a very different behavior. In particular, 
we do not expect to see wave-like behavior in consensus equations, whereas they in our equations 
those play a prominent role. In what follows we will restrict ourselves to coupled second order 
differential equations. 

It is of significant importance to develop a theory that deals with systems where agents may 
interact with few nearby agents. In the case of physical systems with symmetric interactions and 
no damping (such as harmonic crystals), this theory exists and can be found in textbooks [17]. 
It consists of imposing periodic boundary conditions^ and then asserting that the solutions of the 
periodic system behave the same way as in the system with non-trivial boundary conditions, 
except near the boundary. Although we know of no formal proof in the literature that this is 
correct, this method of solution has been used for about a century with great success. Neither is it 
the case that we can rely on previous studies of discretization of a second order partial differential 
equation. Indeed the finite difference method applied to a wave equation with convection will 
give rise (for small enough mesh) to nearly symmetric equations [18,19]. 

In flocks there is no reason for the interactions to be symmetric or undamped, as is the 
case in the study of harmonic crystals. The equations studied here are therefore more general. 
Furthermore in flocks it is desirable to have a two parameter set of equilibria, namely motion with 
constant velocity and constant distance between any two consecutive agents {coherent motion). 
Thus it is necessary to study a more general problem, namely convergence to coherent motion in 
the presence of asymmetry and damping. In this paper we generalize the previously successful 
approach (periodic boundary conditions). In doing this one needs to be aware though that (i) 
assumptions or conjectures needed to solve the old problem must be investigated again as they 
may not be justified anymore, and (ii) new phenomena may arise. For more details see [20,21,9]. 

In the case of linear response theory in solid state physics [17], when a system of symmetrically 
coupled undamped oscillators is perturbed, the signal will typically travel through the entire 
system at constant velocity without damping. In our case, the system is generally either stable 
or unstable. In the former case the perturbation will die out over time, and in the latter, the 
perturbation will blow up exponentially in time. However, even in the stable case perturbations 
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may get very large before dying out. The largest amplitude of a perturbed system that is stable, 
may in fact still grow exponentially in the size of the system. This phenomenon is called flock 
instability. Just like “normal” instability, flock-instability is an undesirable property, since it 
makes large flocks unviable. Flock-instability in arrays of coupled oscillators was illustrated 
in [22], and bears similarity to certain phenomena discovered earlier in fluid mechanics [23, 
24] . Thus the first task is to find criteria to identify those systems that are both stable and flock 
stable. 

We thus need to replace the traditional approach using periodic boundary conditions by 
another that we now outline. For those systems that are stable and flock stable (and only for 
those), we conjecture that for times of length 0{N) (where N is the size of the flock) the solutions 
of the periodic system behave the same way as in the system with non-trivial boundary, except 
near the boundary where additional effects must be taken into consideration. It turns out that 
with those constraints the system with periodic boundary condition behaves like a wave-equation. 
Since the travel time of a wave between the leader (agent 0) and the last agent (numbered N) is 
proportional to we can study the dynamics of the perturbed system for times needed up to 
a finite number of reflections. Due to the asymmetry, wave-packages traveling in the positive M 
direction may have a different signal-velocity than waves traveling in the opposite direction. It 
turns out we can use this effect to achieve either substantial attenuation or magnification of the 
traveling wave at the boundary near agent N. 

In the present work we extend this analysis from nearest neighbor systems done in [20,21] to 
next nearest neighbor (NNN) systems, and in doing that we uncover another new phenomenon. 
We will see that for stable and flock stable systems there are still two signal velocities, but that 
in contrast with nearest neighbor systems it is possible that they have the same sign. This means 
that perturbations can travel (as waves) in only one direction. As a consequence, they cannot be 
reflected. This type of transient has the counter-intuitive characteristic that they travel through 
the system in finite time, after which the system finds itself in (almost) perfect equilibrium. 

The paper is organized as follows. In Section 2 we define the model of interacting agents. The 
main line of reasoning of the method is given in Section 3. The stability conditions are given 
in Section 4. The description of the stable solutions is presented in Section 5. This includes the 
description of the reflectionless waves on the line, which to the best of our knowledge is a new 
result. We include extensive numerical analysis in Section 6 to back up our theory. 


2 The Equations of Motion of the NNN System 

We consider a model of an one-dimensional array of linear damped coupled (up to next nearest 
neighbor) harmonic oscillators on the line. The oscillators or agents are numbered from 0 to A" 
from right to left. We impose that the system is decentralize^ that is: the agents perceive only 
information about other agents that is relative to themselves, in this case relative position and 
relative velocity. See Figure 1 for a sketch of information flow. 

The equations of motion of such a system can be written as: 

2 

^ ^ [QxPxJ (^/c+jf T j ^) T QvPvJ {^k-\-j J'/c)] 5 (f) 

i=-2 

where A is the desired inter-agent distance and pxj {pvj) are position (velocity) parameters. 
The latter are normalized so that px^o = Pv,o = 1 - The normalization factors, px and are often 
called the ‘gains’ in the engineering literature. 
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Px, — 2 Px,2 

Pv,-2 Pv,2 


Fig. 1 Sketch of information flow. Available information about position pxj and velocity weight of 
nearest j = k ^1 and next nearest j = k ^2 agents for fc’th agent. 


The initial conditions we will impose from here on, are as follows. At time t < 0 the agents 
are in equilibrium, Xk = —k A. Then, for t > 0, the leader xq starts moving forward at velocity 
t;o: 


Vt > 0 xo{t) = vot . 

The leader is not influenced by other agents, although other agents (e.g., k = 1 and k = 2) are 
influenced by it. This choice of the leader at the head of the flock is motivated by applications 
to traffic situations (see [20,21]). It is possible to analyze the dynamics with leaders in different 
positions or having more than one leader. We will not pursue this here. 

It is convenient to eliminate the constant A from Equation (1), using the change of coor¬ 
dinates: Zk = Xk ^ k A. In this notation, the equation of motion of the flock in M becomes: 

Definition 2.1 The equations of motion of the NNN system with N > 4 agents, for k G 
{!,••• N}, are: 


2 

’^k — ^ ^ {pxPxJ^k-\-j T PvPvJ^k+j^ • 
i=-2 

This system is suhjeet to the eonstraints 


2 2 

Px,0 — Pv,0 ~ 1 5 ^ ^ Px,j — ^ ^ Pv,j — 0 5 

i=-2 j=-2 


( 2 ) 


( 3 ) 


and to the initial eonditions: 

Zk{0) = 0 , 4(0) = 0 , and ZQ{t) = vot . 

From now on we denote this system by Sn- The eolleetion of the systems {Sn}j^^ 4 ^ will be denoted 
by S. 

Now we use vector notation and write z = {zi, Z2^ Z3, ... zjy)^ together with i = 
(zi, 4, ^ 3 , ••• ‘ Equation (2) may be rewritten as a first order system in 2N dimen¬ 


d 

dt 





sions; 


0 I 


+m. 


( 4 ) 
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where Lx ,Ly G are matrices - the Laplacians - with standard definition 

2 2 

^ ^ 9xPx,j^k-\-j •) ~ ^ ^ 9vPv,j^k-\-j 5 (5) 

i=-2 j=-2 

where F{t) is the “external force” that describes the influence of the leader with trajectory 
zo{t) = vot on the acceleration of its immediate neighbors. It is easy to check that all components 
are zero except the N + 1-st and N + 2-nd components. The exact form of that external force 
depends of the boundary conditions we choose to impose on the system, as we briefly discuss 
now. 

Note that the equations for zi, and zn are subject to non-trivial boundary conditions 

(BC), because there are no agents with numbers — 1, A/" + 1, and N 2. So the equations of 
motion for agents 1, N — 1, and will have to be modified. Here we will use two sets of BC: 
fixed interaction and fixed mass. In the case of fixed interaction BC the central coefficients, px^o 
and of the boundary agents are not equal 1, instead it is the sum of existing interactions. 
On the other hand, in fixed mass BC we change the interactions and keep the central p’s equal 
to 1. The details are given in the Appendix A. 

Coherent motion is defined as: 


yk{t) = aot + bo-kA, (6) 

where ao and bo are arbitrary real constants. It is easily checked that coherent motion is a solution 
to the differential equations given above. Our aims are: 

1. To find out for which values of the parameters trajectories the system is stable: namely, for 
all k, limt^oo \xk{t) ~ yk{t)\ =0 where yk is given in Equation (6). 

2. To find out how fast the stable systems converge to its coherent motion. 

3. To determine what is the size of the transient maxt>o \xN{t) — PAr(t)| in stable systems. 

In the last item we consider only the last (or A^-th) agent to simplify the exposition. As an 
example in Figure 2 we present a sketch of the dynamics expected in the stable system of locally 
coupled oscillators on the line. In the figures we plot the positions relative to the leader, i.e., 
Xk{t) - Vot. 


3 Method 

The analysis of the system of Definition 2.1 is very difficult because the Laplacians given in 
Equation (5) are not simultaneously diagonalizable. In order to overcome that we define a system 
where the communication structure is not a line graph but a circular graph ( see Figure 3). We 
use it (in Section 4) to deduce necessary conditions for stability and flock stability on the line, 
and (in Section 5) to derive expressions for the signal velocities. 

Definition 3.1 The equations of motion of the system with periodic boundary conditions (PBC) 
are: 

2 

’^k — ^ ^ {9xPxJ^k-\-j 4“ 9vPv,j^k-\-j^ • 
i=-2 

This system is subject to the constraints 

2 2 

Px^O — Pv,0 ~ 1 5 ^ ^ Px,j — ^ ^ Pv,j — 0 • 

i=-2 j=-2 
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Finally, instead of boundary conditions for zi, zn-i, cind zn, we set: 

j ^N-\-j — • 

From now on we denote this system by S'ff. The collection of the systems denoted 

by S\ 

The Laplacians L* [with the same definition as in Equation (5)] now become circulant matrices 
and are therefore diagonalizable by the discrete Fourier transform [25]. Let Wm be the m-th 
eigenvector of I/*’s, that is the vector whose j-th component satisfies: 

{W^)j = , 

with (j) = ^'KmjN. We denote the m-th eigenvalues of L* by \x^m and those of L* by With 

a slight abuse of notation we also consider these eigenvalues to be functions Xx{<f) and A^(0) of 
(j) defined above. By using the m-th eigenvector above to calculate L%Wm and L^Wm it is easy to 
show that: 

Lemma 3.1 The A’s are given by 

2 2 

= 5® E = 5a; E cos(j(/)) + sin(j(/))] , 

j=-2 j=0 

2 2 

K{(f>) = 5« E = 5« E [“"’J sin(j</>)] • 

i=-2 j=o 


Stationary orbit position of 

of last agent the leader 



Fig. 2 Dynamics of locally coupled arrays. Sketch of time-dependent dynamics of locally coupled oscillators 
on the line (system Sn) of (a) Type I and (b) Type II (see Section 5 for detailed analysis of these solutions), 
x-axis depicts relative position with respect to the leader. 
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Fig. 3 A circular graph. The numbers indicate the labels of the agents. 


Here we have used the following convenient notation. 

Definition 3.2 Let ax,o = o:v,o = 1 ci'^d = 0- For j > 0 we define: 

^X,j ~ Px,j H“ Px,—j 5 Px,j ~ Px,j Px,—j 5 
^v,j Pv,j “1“ Pv,—j : l^vj Pv,j Pv,—j • 

Note that the sum of the a’s equals 0 by Equation (3). 

Let us now focus on the eigenvectors and eigenvalues of associated with Wm- Denoting 
the eigenvalues by we get: 


( T* /*) ( V 

\ ^x J V ^m,:t / ’ \ ^m,± / 

Thus the evolution of an arbitrary initial condition is given by: 



( 7 ) 

( 8 ) 


where the and bm are determined by the initial condition at t = 0. 

Next, let us evaluate the second row of Equation (7) using that are eigenvectors of L*: 

Lemma 3.2 The eigenvalues of are the roots of the eharaeteristie equation 

- XvicTty - ^, (9) 

where = 2'KmjN. The eigenvalues of S'" are a dense subset of the elosed eurves ^ C 

and U- : ^ C defined by Equation (9). 

Our treatment follows that of [21] where it is conjectured that (for nearest neighbor systems) 
a circular system and a system on the line evolve in a similar manner. The result is that we 
can analyze the circular system and apply the conclusions to the systems on the line. We briefly 
outline how the evolution of the two systems can be compared. 

First we need to remind the reader of the two notions of stability that play a crucial role in 
our analysis. 

Definition 3.3 For given N, the system Sn (S^) is asymptotieally stable if given any initial 
eondition, the trajeetories always eonverge to a eoherent motion and the eonvergenee is exponen¬ 
tial in time. For the systems we eonsider this is equivalent to: Mn (Mf^) has one eigenvalue 
zero with multiplieity 2, and all other eigenvalues have real part (strietly) less than 0. Sn (S^) 
is unstable if at least one eigenvalue has positive real part. 
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Flock stability was introduced in [22]: 

Definition 3.4 The collection S is called flock stable if the Sn are asymptotically stable for all 
N and z/niaxteR|^Ar(^)| grows sub-exponentially in N. 

Note that asymptotic stability is different from flock stability. The former deals with the growth 
of the response of a single system as t tends to infinity while N is held fixed, while the latter 
deals with the growth of the response of a sequence of systems as N tends to infinity. 

Now we mention the main ideas that allow us to compare the evolution of the two systems. 
The first idea is the conjecture that if the system on the circle is asymptotically unstable, then 
the system on the line is either asymptotically unstable or flock unstable. Notice that undamped, 
symmetric systems are all marginally stable, and this aspect does not enter the traditional dis¬ 
cussion in the physics context. This gives us necessary conditions for stability and flock stability 
of the system on the line. 

The second idea involved in this analysis is the principle that, if the system on the line is 
stable and flock stable, then the evolution away from the boundary of the two systems should be 
the same. As we shall see this means that for these systems we obtain wave-like behavior with 
signal velocities determined by the eigenvalues of the system on the circle (see Theorem 5.1). 
This is similar to what is commonly known in solid state physics as periodic boundary conditions 
(see Chapter 21 in [17]), though not exactly the same. The difference is that here we apply 
principle in more generality than is usual in physics, because we are considering systems that 
are not symmetric and not Hamiltonian. 

We know that the reverse of this conjecture is actually false: stability on the circle does not 
imply stability on the line. There are systems that are stable if periodic boundary conditions 
are imposed, but have some eigenvalues with positive real parts when given non-trivial physical 
boundary conditions. In Figure 4 we show a simulation of such a system on the line. The pa¬ 
rameters are given in the Figure. Another example is given in Section 5. It turns out, perhaps 
fortunately, that such counter examples are extremely rare. 

The third and last idea is that the cumbersome physical boundary conditions of Appendix A 
may be replaced by a single “free boundary condition” and a single “fixed boundary condition”. 
This is a great simplification, because the set of possible all physical boundary conditions form 
a 16-parameter set, with no obvious naturally “preferred” boundary condition. However because 
of this last principle, our conclusions will be independent of the physical boundary condition. As 
before, in the traditional physics context, this problem play little or no role, because presumably 
the fixed mass BC is the only possible BC. 

In extending the principle of periodic boundary conditions and adding some new ideas to it, 
we need to be aware that new phenomena may appear (see Section 5.2) and indeed its validity 
is not guaranteed nor is it implied by the validity of the principle in the restricted (symmetric, 
undamped) case (nor indeed by the validity in the general nearest neighbor case). Thus our 
conclusions need to be checked numerically (see Section 6). 


4 Stability Conditions 

We wish to establish conditions that guarantee that the systems Sn on the line is both asymptot¬ 
ically stable and flock stable. Since a direct verification is too hard or even impossible to perform, 
we use the conjectures stated in Section 3. According to those, necessary conditions include the 
stability of the systems 5*, a much simpler problem. 
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a 


N = 200,A = 1,vq= l,gx= -2,g„= -2 
Px = (4/27, -289/432, 1, -253/432,23/216) 
Px = (47/216, -29/108, 1, -79/108, -47/216) 



Relative position 


Fig. 4 Dynamics of unstable system. Dynamics of example system Sn as calculated for N = 200, A = 1, 
VO = l^g^ =g^ = -2, = (4/27, -289/432,1, -253/432, 23/216), pv = (47/216, -29/108,1, -79/108, -47/216) 

and fixed interaction BC. Each color represents the orbit of one of the 200 agents. 


Substituting the expressions for the A’s in Lemma 3.1 into Equation (9), we see that the 
eigenvalues of are the roots of the following equation: 

2 2 

-vgv ^ Pv,3 - 5* E ^ ® do) 

i =-2 j =-2 

Note that when 0 = 0, the characteristic equation becomes = 0. This gives two zero eigen¬ 
values. These trivial eigenvalues are associated with the coherent solutions of the system, Zk = 0 
[see also Equation (6)]. 

Lemma 4.1 The following are neeessary eonditions forS^ not to have eigenvalues with positive 
real part when N is large: 

(i) Px,l + 203,^2 = 0; 

(ii) gv < 0; 

(Hi) ay^i e [-4/3,0], 

(iv) gx(^x,i ^ 0* 


Proof: To prove (i) notice that the roots of characteristic Equation (9) are: 




A^(0) ± y/A^X^p"T4A^^ 


( 11 ) 
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As 0 = 2'KmlN becomes very small, the A’s can be approximated by their first order expansion. 
From Definition 3.2 and Lemma 3.1 we obtain: 

2 2 

\x{4) 0) « tgx<p'^jPx,j, A„((/)-;> 0) « %gv4>'^jPv,j ■ 

i=o i=o 

Substituting these into equation for z^, Equation ( 11 ), we see that for small enough 0, the term 
±A/4Aa,(0) dominates. Since (j) can be either positive or negative, this has four branches meeting 
at the origin at angles of 7 r/ 2 . Two of these branches contain eigenvalues with positive real part 
(for big enough N). Therefore, for N large enough there are 0 such that z^±(0) have negative 
real part unless = 0 - 

For condition (ii) we note that the mean of the two roots of Equation ( 11 ) is equal A^/ 2 . It 
follows that we must require 5R[A^((/))] < 0 for all 0 7 ^ 0 . Since the average ^ 5R[A^(0)]d(/) is 

there is a 0 so that J?[Av(0)] > g^. That of course means that g^ must be non-positive. 

For (hi) note that 5R[A^(0)] < 0. Therefore, ay j cos jcf) > 0. For the NNN system, the 
constraints on the a’s now give 

1 + ay^i cos{(j)) — (1 + cos( 2 (/)) > 0 . 

Since cos(20) = 2 cos^ 0 — 1, the inequality becomes a quadratic inequality in cos(0): 

— (2 + 2a^^i) cos^( 0 ) + ay^i coscj) + 2 + ay^i > 0 , 


which factors as: 


— [(2 + 2a^q) cos[(j)) + 2 + a^q] (cos(0) — 1) ^ 0 . 

By working out three cases, a^q < —1, a^q = —1, and a^q > —1, the conclusion of (iii) may 
be verified. 

Beside 0 = 0, one other case of Equation (9) is easy, namely 0 = tt with the A’s as defined in 
Lemma 3.1 


2 2 

^ ^9v ^ ^ ( 9x ^ ^ ( (^x,j — 0 . 

j=0 j=0 

The roots have non-positive real part if and only if both coefficients are nonnegative. In particular, 
this implies that last term in the above equation is gx {~^y^xj ^ 0. From Definition 3.2 

we know that ^x,j = 1 + Sj=i ^xj = 0, and as a consequence gxOix,i ^ 0, which is 

condition (iv). Similarly, ^^^a^q > 0 but this already follows from conditions (ii) and (iii). □ 

Since we are only interested in the parameter values for which the collection 5* is not unstable, 
we use the above Lemma 4.1 and Definition 3.2 to eliminate a few parameters from our equations. 
This is done by eliminating Px, 2 , <^cc, 2 , and through the substitution: 

yx,2 = —-^Px,l 5 Oix,2 = “(1 + Oix,l) 5 <^v,2 = “(1 + <a^q) , 


which we will use from now on. 
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Proposition 4.1 If the collections'' is stable, the low-frequency expansion o/z^±(0) '1^ given by 

gv{Pv,i + 2 /?„, 2 ) ± f gl{Pv,i + 2/3„,2)^ - 2c/3;(4 + 3q!3;,i) 

,0^ //I I Q ^ ^ QviPv,! ^‘^9xPx,l 

H-r Q?;(4 + oOo, 1) ±— , — 

4 [ ’ \/9l{Pv,i + 2/3^, 2 )^ - 2^^(4 + 3o^,i) _ 

Proof: One can transcribe the first two terms of the corresponding expansion given in [20], or 
one can find the result by substituting power series in 0 in Equation (10) or Equation (11). □ 

This result immediately implies two other necessary criteria for stability. It is unclear whether 
together with the earlier criteria from Lemma 4.1 these also constitute a sufficient set of criteria 
for the stability of 5*. 

Corollary 4.1 The following are necessary conditions for the collection S" to not be unstable: 

('^) 9v{Pv,i + 2/3^, 2 )^ — 2^2,(4 + 3ax,i) > 0, 

dldxi^ + 3(a^q)^(4 + 3ax,i) + ‘^g^gxi^v,! + 2^^^2)(4 + ^oty^i)Px,i + ‘^QxPx,! — 0* 

Proof: If condition (i) does not hold, then one branch of the first order expansion given in 
Proposition 4.1 will have positive real part. Condition (ii) corresponds to setting the argument 
of in Proposition 4.1 as negative. □ 

Remark 4.1 lEe summarize the stability criteria for later use. From Lemma ^.1 and Corol¬ 
lary 4 c we get a list of necessary conditions for system stability. We added condition vii which 
was derived in Corollary 6.1 using Routh-Hurwitz stability criteria (details are given in Ap¬ 
pendix B). 

( 1 ) Px,l + ‘^l^X,2 = 0; 

(ii) gy < 0, 

(hi) ay^i G [-4/3,0], 

(iv) gxC^x,! ^ 0; 

b) gl{Pv,i + Wv,2f - ‘i'gxi^ + Sax,i) > 0 , 

(vi) + 3Q!„^i)^(4 + + 2glgx{l3v,l + 2/?«,2)(4 + 3Q:«,i)/?a:,l + “^gxPx,! — 0, 

(vii) Qx - gl E ]=-2 pI ,3 ^ 0 


IX\ 

V±{ 4 >) = y 


5 Characterization of Solutions 

We assume that we start with an initial condition given as Equation (8). 

Theorem 5.1 Let Kq > 0 fixed. Suppose the collection S" is stable and that the initial condition 
is such that there are a G (0,1) and q > 0 such that NamVn^~^^ and NbmTn^~^^ are bounded, and 
(2 — q)a < 1. Then for large N there are functions /+ and /_ such that the solutions Zj{t) of 
satisfy 

lim sup \zj{t) - vot - f_{j - c_t) - /+(j - c+t)| = 0 . 

The signal velocities c± are given by 


C± = -]^9v{Pv,i + 2 / 3 ^, 2) ± ]^\j9l{Pv,i + 2 ^^, 2)^ - 2^^(4 + 3 oa,q). 
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Sketch of Proof: If 5* is stable then Definition 3.3 and Lemma 3.2 imply that the eigenvalues he 
on curves bounded away from the imaginary axes, except near 0 = 0 where we have an eigenvalue 
0 with multiplicity 2. The low-frequency expansion of iy± (Proposition 4.1 and Corollary 4.1) 
implies that in a neighborhood Jq of 0 = 0 we can write 

^±( 0 ) = '^ 0^±1 + 0^^±2 + • • • : 

where B±i,B ±2 ^ ^ and furthermore B ±2 < 0- For N large enough, none of the eigenmodes 
survive long enough to travel around the system [t of order 0{N)], except those with ^'KvnjN in 
the neighborhood Jq. For these wave-numbers and times scales we may now neglect dissipation. 

We use the initial condition of Equation (8) with bm = 0. Neglecting dissipation, the evolution 
of the j-th component can then be written as 

Zj{t) = . 

m m 

If we write this as /+(j — c+t), we see that c+ = —B-i. Similarly by setting = 0 (instead 
of bm = 0) one shows that c_ = —B_^i. The general case follows by superposition of these two. 
This yields the asymptotic form of z^it). 

To actually prove the remainder indeed tends to zero, one needs the assumption on the decay 
of the am and bm- This part of the argument is given in [20]. □ 


Remark 5.1 The signal velocities c_ and c+ are in units of number of agents per unit time 
(not in distance per unit time). A positive velocity means going from the leader towards the last 
agent. 

Theorem 5.1 states that if 5* is stable, then for large N the systems 5]^ will evolve like a 
wave equation. From the conjectures discussed earlier we conclude that the solutions oi Sn - for 
large N - will behave the same way, except near boundaries. Near the boundaries we apply the 
appropriate boundary conditions ( see below) to get the final solution. This gives linear growth 
of the transients, and that cannot be improved upon. 

If these conditions are not met, in particular if 5* is unstable, then the conjectures tell us 
that S is either unstable of flock unstable. In the first case the coherent motions are unstable 
solutions, and in the second, transients are exponential in N before dying out. 

It turns out that there are several types of wave-like solutions. These depend on the signs 
of the signal velocities c± given in Theorem 5.1 - see the phase diagram presented in Figure 5. 
There are, in principle, three types of wave-like solutions. We study these separately. 

In our analysis below we ignore cases when c± = 0 or c+ = c_. These cases are interesting by 
themselves, but have properties that make them undesirable for situations like traffic and other 
types of flocking. Thus we do not investigate them here. For example when c_ = 0, distances 
between agents don’t tend to the desired distance Z\, but rather to some value that depends on 
the initial conditions. If c+ = c_, which only occurs in Type II solutions, the velocity of the last 
agent is unbounded as N tends to infinity. 


5.1 Type I: Stable, Flockstable, and c_ < 0 < c+ 

When c_ < 0 < c+ the solutions resemble the traditional damped wave reflecting between 
the ends of the flock. The difference in the signal velocities causes the wave to be damped (or 
magnified) when it reflects in agent N. These solutions are called Type I. 
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c_ 



Fig. 5 Phase diagram of signal velocities. Sketch of three type of solutions. Type I and Type II are stable 
with PBC and also on the line. Type III solution is stable only with PBC. 


For these solutions, it can be shown that the orbit of the last agent [see Figure 2(a)] by the 
k-th amplitude Ak^ the period T, and the quotient \ AkJri/^k\ which we refer to as the attenuation 

a. 


Theorem 5.2 Suppose S satisfies the eonditions of Theorem 5.1. If C- < 0 < c+, then for large 
enough N and at time seales t = 0{N), the system has Type I solutions eharaeterized by: 


^/c 


-vqN 

c+ vc+y 



T = 2N 




where A^, ol, and T are defined above, and c± as in Theorem 5.1. 


The proof is essentially that of [20,21] and relies on two insights. The first is that the high 
frequencies die out fast, so that we only need to consider low frequencies (as in the proof of 
Theorem 5.1). The second is that we replace the physical boundary conditions in 5^ by new 
boundary conditions to get the correct reflections at the ends, namely a fixed boundary condition 
at the leader’s end and a free boundary condition at the other end: 


zo{t) = 0 and Z]sf{t) — ZN-i{t) = 0. 


Because for large N only low frequencies survive, these conditions can be replaced by 


zo{t) = 0 and 


dk 




= 0 . 


( 12 ) 
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That means that near the leader, a pulse reflects (with opposite sign), and near the free boundary, 
the traveling pulse is reflected with the same sign and with amplitude multiplied by a factor 
|c_/c+|. The details are written out in [21]. 

In order to get strong damping to minimize transients, we want |c_| < c+. This means 
that in the velocity Laplacian, more emphasis should be placed on the upstream (lower labels) 
information. Such system have asymmetric interactions. 

Corollary 5.1 Suppose S is asymptotically stable and flock stable. Sn has solutions of Type I 
with \c- \ < if: 

(i) gy {fly^i + 2 fly^ 2 ) < 0 ond 

(ii) ^a: (4 + 3o^p) < 0. 

Proof: If 5 is asymptotically stable and flock stable, then all are stable (by our conjectures). 
The conditions c_ < 0 < c+ and |c_| < c+ imply that c_ + c+ > 0. This implies (i). Statement 
(i) together with c_ < 0 implies (ii). □ 

In Figure 6 (parameters are given in the figure) we present typical dynamics of Type I stable 
system Sn- The characteristics predicted from Theorem 5.2 are Ai = 80, o = 0.4, T = 560, and 
from the simulation we measured Ai = 77.2, a = 0.377, T = 568. 


a 


1000 


800 


600 


400 


200 


N = 200, A = 1, 'Co = 1, = -2, 

(-1/2,1/4,1,-3/4,0) 
(-1,3/4,1,-1,1/4) 


= -2 



-200 -150 -100 

Relative position 


Fig. 6 Dynamics of Type I solution. Dynamics of example system as calculated for N = 200, Z\ = 1, 
VQ = 1, Qx = Qv = —2, px = (-0.5,0.25,1,-0.75,0), py = (-1,0.75,1,-1,0.25) and fixed interaction BC. Each 
color represent the orbit of one of the 200 agents. 
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5.2 Type II: Stable, Flockstable, and 0 < c_ < c+ 

When 0 < c_ < c+, that is the signal velocities are both positive, the wave cannot be re¬ 
flected, because it cannot move with negative velocity. We denote these solutions as Type II or 
reflectionless waves. It was proved [21] that such solutions cannot occur with nearest neighbor 
interactions. 

Since both signal velocities are positive, there is no reflection possible at k = N agent. Thus 
the boundary condition at k = N is useless, and we need another boundary condition. We replace 
Equation (12) by the somewhat counter-intuitive condition: 

zo{t) = 0 and zi{t) — zo{t) = 0 . 


As before for large N only low frequencies survive, and so these conditions can be replaced by 

zo{t) = 0 and ^ = 0. (13) 

Thus we have both a free and a fixed boundary condition at the leader’s end. For Type II, the 
orbit of the last agent [see Figure 2(b)] can be characterized by the amplitude A, the first response 
time Ti and the second response time T 2 . 


Theorem 5.3 Suppose S satisfies the conditions of Theorem 5.1. //O < c_ < c+, then for large 
enough N and at time scales t = 0{N), the system has Type II solutions characterized by: 


A = 





(14) 


where A, Ti, and T 2 are as above, and c± as in Theorem 5.1. 

Proof: Ti and T 2 are the (positive) times at which ZN{t) — zo{t) changes velocity. These can 
be deduced from a Proposition whose reasoning is different enough from earlier work, that we 
include a sketch of the proof in Appendix C. A = TiVq is the distance traveled by the leader in 
the time interval [0,Ti). □ 

Corollary 5.2 Suppose S is asymptotically stable and flock stable. Sn has solutions of Type II 
(both velocities positive) if: 

(i) (/3„,i + 2/3„,2) < 0 and 

(ii) 0 < 2g^ (4 + < gl (/3„,i + 2/?„,2)^- 

Proof: Similar to the proof of Corollary 5.1. □ 

In Figure 7 we present typical dynamics of Type II stable system ( parameters given in the figure). 
The characteristics predicted from Theorem 5.3 are A = 43.845, Ti = 43.845, T 2 = 456.16, and 
from the simulation we measured: A = 43.182, Ti = 43.182, T 2 = 453.95. From the figure it 
seems that a start signal traveling with velocity c+ and a stop signal traveling with velocity c_ 
travel from the leader towards the last agent. A striking aspect of this type of solution is that 
very briefly after the second response time, the trajectory of the last agent is (almost) exactly in 
its equilibrium position. Dynamics within such a system can be described as a traveling wave- 
package which does not reflect in the boundary of the system. 
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N = 200, A = 1,1)0= 9x= -2,5„= -2 

p^= (1,-2,1,0,0) 
p„= (-1/2,-1,1,1/2,0) 



Fig. 7 Dynamics of Type II solution. Dynamics of example system as calculated for N = 200, A = 1, 
Vo = 1, gx = gv = —2, px = (1, —2,1, 0, 0), pv = (—0.5, —1,1, 0.5, 0) and fixed interaction BC. Each color represent 
the orbit of one of the 200 agents. 


5.3 Type III: c_ < c+ < 0 

Finally, when c_ < c+ <0, the perturbation which in our set-up starts at the leader, cannot be 
transmitted to the flock, because only negative signal velocities are available. Thus the system 
“finds” another non wave-like solution which has very large amplitudes. The only reason for 
listing this solution in this work at all, is that the system is stable and does have wave-like 
solutions with negative signal velocities. We call these solutions Type III. As with Type II, these 
solutions cannot occur with only nearest neighbor interactions. 

Corollary 5.3 Suppose 5* is asymptotieally stable. Sn has solutions of Type III (both veloeities 
negative) if: 

(i) gy {^y^i -h 2(3y^2) > 0 and 

(ii) 2^3, (4 + 3Q(a,,i) > 0. 

Proof: Similar to the proof of Corollary 5.1. □ 

Within such a setup, on short time scales, the leader simply starts and other agents do not 
follow him. On time-scales larger than 0{N)^ other phenomena may take place. Thus amplitudes 
will grow faster than 0{N), and the system is flock unstable or even asymptotically unstable. 
However, due to the complicated nature of the stability conditions, we do not have a proof of 
this. In Figure 8 we present a simulation of such a system ( parameters given in the figure). 
Notice that the amplitudes do not grow linearly with system size. 
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A = 1, '?;o = 1, = -2 , = -2 

p, = (-2,-15/4,1,-21/4,5/2) 

A = 50 = (-1,4,1,-5,1) A = 100 



Fig. 8 Dynamics of Type III solution. Dynamics of example system Sn as calculated for A = 1, i;o = 1, 
9x = Qv = —2, px = (—2,—15/4,1,—21/4, 5/2), pv = ( —1, 4,1, —5,1), fixed interaction BC N = 50 (left panel) 
and N = 100 (right panel). Each color represent the orbit of one of the 50 and 100 agents, respectively. 


6 Numerical Tests 

As we saw in Section 5, measured values of certain characteristics presented for N = 200 differ 
slightly from the predicted ones, given by Theorem 5.2 and Theorem 5.3. This is expected, since 
our predictions are valid for A" ^ oo. In order to test our conclusions we performed extensive 
numerical calculations. We outline our procedure. 

First we fixed Qx = 9v = and defined a set P of about 8.6 * 10^ ten-tuples 
{9x, 9v, Px-2, Px-I, Px,i, Px,2, Pv-2, Pv-I, Pv,i, Pv,i, Pv,2) Satisfying Equation (3). We call these 
ten-tuples configurations. From this set of configurations we then selected the set Pc that satisfy 
all the criteria in Remark 4.1. For Type I solutions we impose an additional constraint, namely: 
|c_| < c+ (see Corollary 5.1). Next, from the same ten-tuples of configurations we created the set 
Ps^N that satisfy Definition 3.3 for given N. It turns out that for N large enough these sets were 
identical: Pc = Ps,n (ia our case we had to go up to A = 800 for a few systems). This strongly 
suggests that indeed the criteria in Remark 4.1 (plus Corollary 5.1) are a very good indicator of 
asymptotic stability of the system on the circle. 

In order to decrease computation time for large A, we imposed a further constraint on Ps that 
selected 500 configurations of Type I and 500 configurations of Type II. The constraints were for 
the period, namely T < 0(10A) (Type I), and for the second response time, namely T 2 < 0(10A) 
(Type II). We ran each of these configurations for A G {25 • that is: for A varying from 

25 to roughly 52, 000. We measure the characteristics directly from numerical simulations and 
compare them with predictions of Theorem 5.2 and Theorem 5.3. For the numerical work we 
used the ordinary differential equation solver of the Boost library [26,27] in a parallel computing 
environment. 

In Figure 9 we present the relative error=\measured — predicted\/\predicted\ of the following 
quantities: for Type I solutions, the first amplitude Ai, the period T, and the attenuation a, and 
for Type II solutions, the amplitude A and the first and second response times Ti and T 2 . We 
plot both the error average (for 500 measurements/configurations) and the worst (largest) error. 
We repeated this experiment for two different types of physical boundary conditions (denoted by 
fixed interaction and fixed mass (Appendix A) to make sure that these did not make a difference. 
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System size N System size N 


Fig. 9 Relative error size scaling. Size N dependence of average and maximal relative error of Type I and 
Type II solutions for two different boundary condition as calculated for N = 25,..., 51200 agents. Notice that 
the plot has log-log scale, therefore slope corresponds to the power of the decay. 


As is clearly visible in Figure 9, the relative errors decrease as N grows, as is predicted by 
the theory. Our numerical analysis is consistent with the statement that - with the exception 
of period T for Type I orbit - the error decreases as 0{1/\/N). The error in the period T (for 
Type I) appears to decrease as 0{1/N). 


Acknowledgements We acknowledge support by the European Union’s Seventh Framework Program FP7- 
REGPOT-2012-2013-1 under grant agreement no. 316165. 


Appendix A: The Physical Boundary Conditions 

We will introduce two sets of boundary conditions for Sn (the system on the line). We per¬ 
formed numerics with both types of boundary conditions in order to support our conclusion that 
for stable and flock stable systems the trajectories are independent of the physical boundary 
conditions. 

Let tSw be the system in Definition 2.1. In decentralized systems the row sum of the Laplacians 
equals 0, that is: Pv,j = 0- This implies that for the system the equations of 

agents /c = 1, A — 1, and N have to be modified. In the case of fixed interaction BC the masses, 
px,o and py^o^ of the agent are not equal 1, instead it is the sum of existing interactions. On the 
other hand, in fixed mass BC we change the interactions of existing agents and keep the central 
px,o and py^o equal to 1. Here are the details: 
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Definition 6.1 (i) fixed interaction BC: 

~ {QxPx^ — I^O “ 1 “ 9vPv, — i^o) 

~ [px {Px,-1 + PxP + px,2) + pv {Pv,-1 + pv,l + pv,2) Zi] 

2 

+ {pxPxJ^l+j + PvPvjZi^j) , 

i=i 

-1 

'^N—1 — ^ ^ {pxPxJ “ 1 “ 9vPv,j 

J = -2 

~ [Px {px,-2 + Px,-1 + Px,l) ZN-1 + Pv {Pv,-2 + Pv,-1 + Pv,l) ^N-l] 

+ {pxPxP^N PvPv,lZ]Sf) , 

-1 

’^N — ^ ^ {pxPx,j^N-\-j H“ PvPv,j^N-\-j) 

i=-2 

“ [px {px,-2 + Px,-l) + Pv {Pv,-2 + Pv,-l) ^n] • 

(a) fixed mass BC: 

= [px {px,-2 + Px,-l) ^0 + Pv {Pv,-2 + Pv,-l) ^o] 

2 

+ {pxpxj^l+j + PvPvjZi^j) , 
j=0 
0 

’^N—1 — ^ ^ {pxPxJPvPv,j 

i=-2 

+ [px {Px,l + px,2) Zn + Pv {pv,l + Pv,2) Zn] 5 

0 

'^N — ^ ^ (yPxPxJ^N-\-j H“ PvPv,j^N-\-j^ 

i=-2 

+ [px {Px,l + Px,2) Zn + Pv {pvp + px,2) Zn] • 


If we use vector notation, the influence of leader on agents 1 and 2 is formulated as an external 
force. Write z = (zi, ^ 2 , ... and i = (ii, Z 2 ^ is, ... iw )^- The equation of motion 

can be rewritten as a first order system in 2N dimensions: 


_d 

dt 






Those terms in the full equation of motion that contain zq or io are written as external force. 
So all components of the external force F are zero except the N + 1-st and N + 2-nd. These two 
components are given by: 

f PxPx, — lZo T pvPv, — lZo 

\PxPx, — 2Zo T pvPv, — 2Zo 
if we impose fixed interactions BC, and 

f Px {Px,-2 + Px,-l) Zq + pv {pv,-2 + Pv,-l) io 
\ PxPx,-2Zo F pvPv,-2Zo 




if we impose fixed mass BC. 
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Appendix B: The Routh-Hurwitz Stability Criteria 

The Routh-Hurwitz criterion is a standard strategy to derive a concise set of conditions that 
is equivalent to the condition that all roots of a given polynomial have negative real parts. In 
various systems similar to the ones discussed here, this criterion gives good results [20,9]. In our 
current case the resulting equations are too complicated to give us much information and we 
only get one more necessary condition for stability that we can use, namely Corollary 6.1. Our 
discussion is based on Chapter 15, Sections 6, 8, and 13 of Ref. [28], where more details can be 
found. 

Theorem 6.1 (Routh-Hurwitz) Assume that the determinants given below are nonzero. Given 
a real polynomial R = + a 2 X^ + aix + oq, all roots of R have negative real part if and 

only if all determinants of the upper-left submatriees (the leading prineipal minors) of: 


/as 

ai 

0 


1 

a2 

ao 

0 

0 


ai 

0 

\ 0 

1 

Ci2 

ao/ 


are positive. That is: as > 0, uq > 0, asa 2 — ai > 0 , and a^a 2 ai — — a^ > 0. 

An equivalent but less well-known set of conditions is given in the following: 

Theorem 6.2 (Lienard-Chipart) Assume that the determinants in Theorem 6.1 are nonzero. 
Given a real polynomial R = x^ ^ a^x^ + a 2 x‘^ + aix + uq; all roots of R have negative real part 
if and only if as > 0, 02 > 0 , oq > 0, and asa 2 ai — a‘^ao — > 0. 

The characteristic polynomial Q of Equation (9) can be turned into a polynomial with real 
coefficients 


R = QQ* = u^- 23?(A„)z/3 + [|A„|2 - 23?(A^)] 

+ 2 [5R(Aa;)3?(A„) + S(A2;)S(A„)] u + \\x'^ ; 

by taking its product with its complex conjugate. Clearly, all roots of Q have negative real part 
if and only if the same is true for R. Notice that in each of the two criteria, one of the equations 
is trivially satisfied, namely ao > 0 (where we are assuming nondegeneracy). Therefore, in the 
Routh-Hurwitz case three equations are obtained. The first two are: 

5R(A^) <0, (15) 

5R(A„) [|A„|2 - 25R(A,)] - [SR(A,)5R(A„) + 9(A,)9(A„)] > 0. (16) 

The third inequality we do not utilize, since it is extremely complicated containing fifth order 
terms. We are left with the above two, which are now necessary conditions for all roots to have 
negative real part. 

Similarly, the Lienard-Chipart stability criterion also gives two necessary conditions for all 
roots to have negative real part: 


5R(A,) <0, (17) 

2SR(A^)-|A„|2<0. (18) 

The third inequality is the same as before and will not be utilized. Since the second inequality 
of the Lienard-Chipart conditions seems less complicated than the corresponding one of the 
Routh-Hurwitz conditions, we will continue with the former. 
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Substituting the expressions for the A’s (Lemma 3.1) we get: 


9v 


2 

^ cos{j(p) 

j=0 



2 


/ 

2 

2 

2 

9x 

y] cos{j(j)) 

-gl < 


y] cos(j» 

- 

y] Pvj 


J=0 



J=0 


J=0 


< 0 , 

< 0 . 


These are complicated relations therefore we will use the equivalent relations averaged over (j). 
The first of these equations was already used in Lemma 4.1. After some calculations we can work 
out the average over 0 of the second relation. This gives the final necessary condition for all roots 
to have negative real part. 


Corollary 6.1 The following is a necessary condition for the collection 5* to not he b unstable: 

2 

9x-gl Tj ® • 

i=-2 


Appendix C: Analysis of Type II Trajectories 


Proposition 6.1 Let Kq > 0 fixed. Suppose that 5* is stable and that 0 < c_ < c+. Suppose 
further that there are o G (0,1) and q > 0 such that and are hounded, and 

(2 — q)a < 1. Then 

lim sup \zN{t)— Z]sf{t)\ = 0 . 

^^^te[0,KoN] 

where Z]sf{t) is given by 


ZN{t) = < 


—t 


c+ 


C-I-—c_ 

0 


)(‘-S) 


t G 
t G 
t G 



The signal velocities are as in Theorem 5.1. 

Sketch of Proof: We consider the equations of motion for the acceleration of agent k. These 
are given by the second derivative with respect to time of Definition 2.1. In those equations, the 
only expression that depends on time is the initial condition of leader. So nothing changes, except 
that now fo{t) = 6{t) (for > 0), where 6 is the Dirac function. We replace the Dirac function by 
a smooth pulse p{t) that enables us to satisfy the decay constraint on the decay of Om and bm 
but with the condition that j p{s) ds = 1. So now we obtain: 

^o{t)=p{t) (19) 

Theorem 5.1 now implies that in 5* we have 


ik{t) = f+{k - c+t) + f-{k- c-t) 


( 20 ) 
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By the periodic boundary conditions conjectures, we see that away from the boundaries the 
behavior of S and 5* is the same. So we have the above relation from t = 0 until the signal runs 
into the boundary at N. 

Setting k = 0 in the last equation and comparing with Equation (19) gives 


p(.t) = f+i-c+t) + f-i-c-t) ■ 


( 21 ) 


The second part of Equation (13) then gives: 


f+{-c+t) + = 0. 


Integrate with respect to t to get 

-f+i-c+t) - —f-i-c-t) = 0 ^ /-(«) = -—/+ (— s) . 

c+ c_ c+ Vc- / 

Substitute this into Equation (21): 

p{t) = /+(s) = — — — p ( —) . 

c+ C+-C- \c+/ 

Now use both of the last equations to eliminate /_ and /+ from Equation (20): 


aw 


c+ — c_ 


P 



C- 


c+ — C- 


p 



Now set k = N and integrate twice with respect to time and add a Galilean transformation 
so that for small positive t we get Z]v(t) = —t. With some rewriting this gives the final result. 
(As before, to actually prove the remainder indeed tends to zero, one needs the assumption on 
the decay of the am and hm- This part of the argument is given in [20].) □ 
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